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ABSTRACT 

Investigations  Into  the  behavior  of  the  gas  flow  behind 
spherical  or  cylindrical  blasts  have  shown  that  secondary  shocks 
arise  within  the  original  detonation  gases.   The  secondary  shock, 
at  first  weak J  is  carried  outward  with  the  expanding  gases. 
Subsequently  it  strengthens  and  bends  back  toward  the  origin^, 
arriving  there  with  high  intensity. 

ES7  using  some  recently  developed  techniques  in  shock 
dynamics  a  theory  is  developed  by  which  the  motion  of  the  main 
shock  wave,  as  well  as  the  formation  and  subsequent  motion  of 
the  secondary  shock,  are  given  by  explicit  formulae.   A  formulation 
for  giving  the  location  of  the  contact  surface  between  the  deto- 
nation gases  and  the  outside  atmosphere  has  also  been  obtained. 
The  results  of  a  specific  problem,  which  has  been  solved  by 
numerically  integrating  the  total  equations  of  motion,  and  also 
checked  experimentally,  are  compared  with  the  results  of  the 
present  theory. 
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NYO-9552 
A  SIMPLIFIED  DESCRIPTION  OP  SPHERICAL  AND  CYLINDRICAL  BLAST  WAVES 

1.  Introduction. 

The  purpose  of  the  present  analysis  Is  to  demonstrate 
the  possibility  of  obtaining  a  satisfactory  description  of 
the  initial  phases  of  spherical  and  cylindrical  blasts,  without 
having  to  go  through  a  nvimerical  Integration  of  the  Eulerlan  or 
Lagrangian  flow  equations.   To  accomplish  this,  some  recent 
developments  in  shock  dynamics  are  used  and  extended  where 
necessary.   In  addition,  a  technique  has  been  devised  by  which 
it  is  possible  to  give  the  motion  of  a  contact  front  in  terms  of 
flow  properties  carried  to  the  front  by  the  characteristics  meeting 
it.   It  is  hoped  that  some  of  these  results  may  be  applicable  to 
the  ntunerical  integration  of  general  flow  equations,  whose 
solutions  contain  moving  boundaries  such  as  shocks  and  contact 
fronts. 

Recent  theoretical  and  experimental  studies  of  spherical 
and  cylindrical  blast  waves  [Ref .1,2,3,7,8]  enable  us  to  give 
the  following  qualitative  description  of  these  processes.   Assume 
at  time  t  =  0  a  gas  sphere  (or  cylinder)  of  radius  x  =  Xq 
under  high  Internal  pressure,   p  =  Ph,   to  be  in  a  still  air 
atmosphere  at  pressure  p  =  p^,  p^  »  p^.   The  names   "gas" 
and  "air,"  respectively,  will  henceforth  be  used  to  describe 
the  fluids  Initially  within  and  outside  of  this  sphere.   For 
t  >  0  an  equalization  (explosion)  process  takes  place.   At  the 
initial  point  t  =  0,  x  =  x^  the  region  x  >  0,  t  >  0  can  be 
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separated  into  five  domains  (Plg.l):  (O) ,  undisturbed  alrj  (l), 
air  which  has  been  overtaken  by  the  main  blast  wave;  (2),  gas 
under  the  influence  of  the  air;  (3),  expanding  gas  uninfluenced 
by  the  outside  air;  (4),  gas  not  yet  disturbed  by  the  expansion. 
The  gas  In  region  (2)  and  the  air  In  region  (1)  are  separated  by 
a  contact  front.   Also,  a  secondary  shock  develops  and  separates 
regions  (5)  and  (2). 

The  Interesting  phenomenon  of  a  secondary  shock  appears 
only  In  spherical  (and  cylindrical)  flows  and  does  not  arise  In 
one  dimensional  shock  tube  studies.   For  the  latter  the  main  shock 
and  the  expansion  come  Into  an  Instantaneous  equilibrium,  being 
separated  by  a  region  of  uniform  pressure  and  velocity.   Physical 
reasoning  given  In  the  literature  [e.g.,Ref.l]  for  the  secondary 
shock  formation  Is  that  the  high  pressure  gas  upon  passing  through 
a  spherical  rarefaction  wave  must  expand  to  lower  pressures  than 
those  reached  through  an  equivalent  one-dlmenslonal  expansion. 
This  is  due  to  the  Increase  of  volume.   Because  of  this  "over- 
expansion",  the  pressures  at  the  tall  of  the  rarefaction  wave 
are  lower  than  the  pressures  transmitted  back  by  the  main  shock 
and  a  compression,  or  second  shock,  must  be  Inserted  to  connect 
the  two  phases.   Ejy  more  mathematical  reasoning  the  situation  Is 
further  clarified:  the  centered  expansion  is  accomplished  through 
the  negative  characteristics  (slope  ^   =  u  -  a  where  u  is 
particle  velocity  and  a  is  sound  speed).   The  characteristics 
at  the  head  of  the  wave  point  in  the  decreasing  x  direction 
(subsonic  flow).   They  then  fan  around  to  the  increasing  x  direction 
as  the  velocity  Increases.   Negative  characteristics  are  also 
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reflected  from  the  main  shock,  but  as  this  shock  expands  It 
becomes  weaker  causing  these  characteristics  to  incline  more 
and  more  toward  the  decreasing  x  direction.   Hence, 
characteristics  of  the  same  family,  but  arising  from  two 
different  sources,  will  tend  to  intersect.   This  is  prevented 
by  the  secondary  shock.   A  similar  situation  arises  [Ref .6]  in 
one-dimensional  flows  when  a  shock  moving  down  a  uniform  tube  meets 
an  area  change.   The  shock  strength  is  attenuated  and  characteristics 
reflected  from  the  shock  will  be  bent.   They  will  meet  the  charac- 
teristics of  the  same  family  which  originated  in  the  vinlform 
section  of  the  tube  unless  a  secondary  shock  is  inserted. 

This  reasoning  with  characteristics  leads  directly  to  a 
method  for  explicitly  determining  the  path  of  the  secondary  shock. 
Near  the  initial  point,  the  intersecting  characteristics  will  have 
nearly  the  same  slope.   It  Is  easily  shown  [Ref .5>p.l59]  that  the 
slope  of  a  weak  shock,  at  each  point,  is  equal  to  the  average  of 
the  slopes  of  the  Incoming  characteristics  at  that  point »   Using 
this  fact,  Whltham  [Ref. 9]  developed  a  technique  for  obtaining 
the  differential  equation  for  the  shock  path.   This  technique  is 
used  with  the  present  problem  to  get  the  initial  motion  of  the 
secondary  shock. 

As  the  secondary  shock  moves  outward  It  strengthens,  and 
the  above  "weak  shock"  theory  is  no  longer  applicable.   In  order 
to  follow  further  shock  development  a  relation  obtained  by 
Chlsnell  [Ref.^]  and  amplified  by  Whltham  [Ref. 10]  is  Introduced. 
Chlsnell  considered  a  shock  wave  moving  In  a  channel  with  varying 
cross  section  area  and  derived  a  closed  form  expression  which 
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gives  the  relation  between  the  channel  area  and  shock  strength. 
By   a  suitable  choice  of  channel  shape  the  analysis  may  be  used 
to  describe  a  spherical  shock.   Whitham  considerably  simplified 
and  extended  Chisnell's  analysis.   He  showed  that  the  "shock- 
area"  relation  can  be  obtained  by  application  of  the  following 
rule:  apply  the  differential  relation  which  is  to  be  satisfied 
by  the  flow  quantities  along  a  characteristic  to  the  flow 
quantities  Just  behind  the  shock.   Using  the  shock  relations 
which  express  these  values  in  terms  of  the  shock  strength, 
another  differential  equation  is  obtained  relating  shock  strength 
change  and  area  change.   For  the  present  analysis  the  channel 
area  is  identified  with  the  shock  surface  area  and,  hence, 
the  shock  radius.   We  obtain,  finally,  a  differential  relation 
between  shock  radius  and  strength.   If  the  flow  ahead  of  the 
shock  is  uniform,  this  relation  can  be  integrated,  giving 
Chisnell's  result.   The  inward  facing  secondary  shock,  however, 
has  ahead  of  it  the  non-uniform  expajidlng  gas.   Hence,  the 
shock  motion  depends  not  only  on  its  radius  but  also  on  the 
flow  conditions  immediately  ahead  of  it.   Whitham's  derivation 
of  the  "shock-area"  relation  is  extended  to  include  this  situation. 

The  secondary  shock  is  therefore  obtained  in  two  sections 
which  take  into  account  its  changing  nature  from  very  weak. 
Initially,  to  quite  strong  when  it  finally  Implodes  at  the 
origin.   When  it  is  weak,  the  average  slope  technique  will  be 
used  to  describe  its  motion  until  the  point  where  it  either 
starts  bending  back  toward  the  origin,  or  where  the  shock  strength 
becomes  too  large.   Prom  this  point  on,  the  extended  "shock-area" 
rule  is  used. 
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The  determination  of  the  main  shock  wave  separating  regions 
(0)  and  (l)  Is  much  simpler  than  that  of  the  secondary  shock. 
This  Is  because  the  flow  ahead  of  the  main  shock  In  region  (O) 
Is  uniform,  and  direct  application  of  the  "shock-area"  relation, 
without  extensions,  suffices.   As  a  consequence  of  the  form  of 
the  "shock-area"  relation  the  x  and  t  coordinates  of  the 
main  shock  are  determined  parametrlcally  as  fxonctlons  of  shock 
Mach  number.   Similarly,  since  the  contact  front  and  secondary 
shock  are  determined  by  extrapolation  of  the  known  flow 
conditions  at  the  main  shock,  their  coordinates  are  also 
determined  as  parametric  fvmctlons  of  the  main  shock  Mach 
number.   Therefore,  starting  at  the  Initial  point  with  the 
known  main  shock  Mach  number,  small  Incremental  steps  In  this 
Vlach   number  are  made.   Corresponding  to  each  Increment,  a  point 
on  the  main  shock,  contact  front,  and  secondary  shock  are  determined. 

Throughout  the  present  analysis,  the  effects  of  entropy  changes 
are  neglected.   This  does  not  affect  the  calculations  in  the  ex- 
pansion region  {J>)    since  the  expansion  Is  Isentroplc.   As  the 
secondary  shock  starts.  It  Is  weak,  and  entropy  changes  In  region 
(2)  are  small  near  the  origin.   In  the  application  of  the  shock- 
area  rule,  it  Is  seen,  using  Whltham's  derivation,  that  the 
characteristic  differential  relation  used  does  not  contain  entropy 
terms  explicitly.   Also,  the  flow  quantities  appearing  there  are 
evaluated  Immediately  behind  the  shock  by  use  of  the  Hugoniot 
shock  relations.   Hence,  if  the  flow  ahead  of  the  shock  is 
isentroplc,  entropy  considerations  will  not  affect  the  shock 
motlonj  this  is  the  case  for  both  the  main  and  secondary  shock, 
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The  contact  front,  being  a  particle  path,  carries  a  constant 
entropy,  or,  as  in  the  present  problem,  the  change  in  entropy 
across  the  front  is  constant.   It  is  seen,  then,  that  the 
only  place  where  neglect  of  entropy  changes  could  affect  the 
calculations  is  in  the  extrapolation  of  flow  properties  from 
the  main  shock  through  the  small  region  (l)  to  the  contact  front. 

In  Section  2  an  approximate  solution  to  the  flow  in  the 
expansion  region  (;5)  is  determined.   Since  the  flow  near  the 
initial  point  can  be  accurately  described  by  one  dimensional 
shock  tube  theory,  the  full  solution  is  approximated  by  perturbing 
about  the  known  one  dimensional  solution.   Equations  describing 
the  main  shock  and  negative  characteristics  in  region  (l)  are 
derived  in  Section  3«   In  Section  H-   the  contact  front  and 
negative  characteristics  in  region  (2)  are  determined  by  the 
utilization  of  known  fluid  dynamical  properties  of  contact 
discontinuities.   The  secondary  shock  equations  are  obtained 
in  Section  5=   Here,  both  the  "weak-shock"  theory  and  the 
extended  "shock-area"  rule  are  presented.   In  Section  6  a 
specific  problem  is  treated,  and  results  are  compared  with  other 
solutions  to  this  problem. 

The  present  analysis  will  describe  the  motion  of  the 
secondary  shock  as  it  initially  moves  outward  and  until  it 
finally  implodes  at  the  origin.   Beyond  that  time,  its  motion 
cannot  be  given  by  this  theory.   Similarly,  this  time  will 
roughly  delineate  the  region  of  validity  for  the  main  shock  and 
contact  front  approximations. 
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2.  Expansion  Region. 

Immediately  after  detonation,  the  high  pressure  gas  in 
region  (k)    of  Figure  1  Is  penetrated  by  a  centered  expansion 
fan  and  moves  rapidly  outward.   This  expansion,  region  (5), 
Is  Isentroplc  and  satisfies  the  flow  equations  (1)  and  (2), 
[Ref.5]*  written  In  characteristic  form 

(1)  (^  ^  *  I  u,)  +  (u+a)  (^  a^  +  I  uj  +  fHa  .  0 

(2)  (ft  S  -  I  "t)-("-*)  (ft  ^x  -  I  "x)  +  55^  =  0 

n  =  1  for  cylindrlcally  symmetric  flows, 
n  =  2  for  spherically  symmetric  flows, 
and  Y  Is  the  ratio  of  specific  heats. 

The  flow  variables  a  and  u,   sound  speed  and  particle  velocity, 
as  well  as  the  Independent  variables  x  and  t,   radial  distance 
and  time,  are  all  assumed  to  have  been  made  dlmenslonless  as  follows; 

a  =  VaQ,  u  =  y^Q,     X  =  Vxq,   t  =  ^^0/xq. 

Barred  quantities  are  true  physical  variables,   a^  Is  the  sound 
speed  outside  the  blast  wave  In  region  (O)  of  Figure  1  and  Xq 
Is  the  Initial  blast  radius. 

An  approximate  solution  to  equations  (l)  and  (2)  Is  obtained 
by  perturbing  about  a  one  dimensional  centered  expansion  wave, 
whose  solution  can  be  given  explicitly  [Ref  .5j.p«104] .   This 
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approximation  holds  exactly  as  the  Initial  point  x  =  1,  t  =  0 
Is  approached.  To  obtain  this  approximation  we  first  make  the 
following  definitions s 


u  =  u,  +  Up 


a  =  a,  +  ag 


O)        <   ^  =  :ct  ^1  +  P  ^1  *  ^ 


a,  -  TT  n^ 


^7^1  »i  -^  2  "1  *  "   7-1  1   2  --1 
-1^  &^  +  ^  xx^   ,  S  =  :^  a2  -  ^  u^. 


The  variables  with  subscript  1  are  solutions  to  the  one  dimensional 
centered  expansion  problem  (equations  (l)  and  (2)  with  n  =  O) : 


u,  =  2nr  +  (l-^i)  ^ 


(4)     ;  a^  =  tx(2r  -  iS^)  j  ^,  =  1^1 


x-1 


r  =  constant,   s  =  (l-2^.)r  -  (l-pi)— r 


Substituting  (5)  and  (k)   Into  (l),  retaining  lowest  order  terms, 
we  obtain  an  equation  for  R% 


nu,  a. 


or 

(5)       li  _       dx  -2xdR 


T 


4pjt+(l-2M.)(x-l)   ntu^a^^ 


The  first  relation  In  (5)  gives  for  the  curvilinear  positive 
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characteristics 

(6)       iS^  .  2r  =  -  Kt■^^^,   or  a^t^^^  =  Kn 

where  K  is  a  constant  for  each  characteristic. 

Using  (6)  to  eliminate  x,  we  obtain  the  differential 
equation  for  R  along  the  positive  characteristics. 


(7) 


dR  ^  -nM.K[2r-(l-pL)Kt"^'^]t"^^^ 
^      2[l+t(2r-Kt"^^^)] 


It  does  not  seem  possible  to  integrate  this  equation  for  arbitrary 
ji.   However,  for  certain  specific  values  a  closed  form  integration 

is  possible.   The  integration  for  M-  =  k  ^^^     IT  *  corresponding 

7       5 
to  y  =  -k     ^^^    "^  »      *^^^  ^®  carried  out,  but  with  the  present 

approximate  theory  the  added  complexity  of  the  exact  integral 

does  not  seem  Justified.   The  expression  on  the  right  hand  side 

of  (7)  is  therefore  simplified  by  assuming  the  bracketed  term 

in  the  denominator  equal  to  1.   This  implies  we  remain  close 

to  the  line  x  =  1,   which  is  roughly  correct  except  for  very 

high  intensity  blast  waves.   With  this  assximption  we  have 


dR  ^  .  n^  [2r-(l-n)Kt-2M.]t-2li. 


This  is  integrated  to  give 


R  =  -  2^[.^  t-^^  -  ^^   t-"-^]  +  f  (a^tS'') 
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Eliminating  K  by  means  of   (6) 

The  arbitrary  function  f(aj^t  ^) ,   constant  along  each 
positive  characteristic,  and  the  constant  r  are  determined  by 
continuously  connecting  -ij  a  +  |-  u  across  the  boundary 
characteristic  between  regions  (5)  and  (4).   If  (4)  were  a 
region  of  uniform  flow  with  u^^  and  ai^  constant,  the  boundary 
characteristic  would  be 

(9)  x-1  =  (u2^-a2^)to 

For  this  case  we  would  have 


r  =  ;pX  ^4  +  I  ^4 


(10) 


„.   .2pL>     tn  /  Ix^  r   1-H   a  2    2r  ^  , 
f  (a^t  n    -  -  -^   (a-)    V(i4n)  ^4  -  T:^  ^4^ 


The  second  equation  In  (lO)  Is  obtained  by  requiring  R  to 
vanish  at  the  boundary  between  regions  (5)  and  (4).   We  shall 
assume  for  deflnlteness  in  later  calculations  (and  simplicity) 
that  region  (4)  is  uniform  and  that  equations  (8)  and  (lO)  hold 
In  region  (3).   It  is  clear  that  by  properly  choosing  the 
arbitrary  characteristic  function  f(a,  t^)   appearing  in 
equation  (8),  a  more  general  boundary  condition  could  be  satisfied. 
Combining  (8)  and  (10)  we  have,  for  region  (3) 

-  1?  - 


a,  "^^ 


with  ,  , 


«(-i)  ff  ^Ti^t^i'-(lt)\'i  -  Tiii[-i-(lt'  %^ 


^1  =l^(2''  -T^)- 


Substituting  the  relations  (5)  and  (4)  Into  the  negative 
characteristic  equation  (2),  and  retaining  lowest  order  terms 
an  equation  for  S  Is  obtained: 


nu,  a,  t 
(12)      tS^  +  (x-l)S^  +  [(2n-l)R+S]  +   g/   =  0, 


where  R  =  tH(aj^),   and  H  Is  defined  In  (11).   Writing  (12) 
In  characteristic  form. 


/,,N      dt  _  dx  _  dS 

^■^^'  ~r  ~  x^ nu  a,  t 

-[(2n-l)R+S] ^3^ 


The  first  of  equations  (13)  give  for  the  characteristics  the 
lines  x-1  =  Lt,   where  L  Is  constant  along  each  characteristic, 
Along  the  characteristics,  (15)  also  gives 

,  nu,  a,  t 

d  /i.o\      /-,      r^     \r,  11 


(It)      ^(tS)  =  {l-£p.)R  -  ^^^) 


nu,a,  t 


=  (l-2M.)tH(a^)  -  ^-^^) 
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since  u   and  a,   are  constEint  along  the  characteristic 
lines,  (14)  can  be  Integrated  Immediately  to  give 

s  .  (1-a.)  I  H(a,)  -  ^  U  -  M^) 

Letting  L  =  (x-l)/t  and  using  equation  (4)  we  have 

(15)  ;^  a^  -  I  u^  =»  (l-2n)r  -  (l-^i)  ^   +  ^^^   tH(a^) 

nu,  a^  t 

^^^^r,   [x-1  -  log  x] 

2(x-l)2 

For  determination  of  the  secondary  shock  It  Is  necessary  to 
have  the  equation  for  the  negative  characteristics.   The  lowest 
order  approximation  for  this  equation  Is  x  =  1  +  Lt.   However, 
by  using  the  above  results,  a  better  approximation  can  be  obtained. 
To  get  this  we  first  obtain  u  -  a  by  taking  a  proper  linear 
combination  of  (11)  and  (I5) 

(16)  u  -  a  =  iy^  +  |l|H  tH(a, )  +  "^^^ ^  [x-1  -  log  x] 

^    ^    ^    ^-^^  ^  2(l-|x)(x-l)2 

Along  the  negative  characteristics  we  have 


dx 

^  =  u  -  a. 


To  facilitate  Integration  of  the  above,  the  right  side  of  (I6) 
Is  simplified  by  assuming  x  =  1  +  Lt,   and  u^,  a^,   and  L 
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are  constant  along  the  negative  characteristics,  giving 


nu,  a. 


dx       T  ^  l-2tJL  ^„/o    ^   J.       11        n        log(l+Lt)i 


nu,a. 


V  -  1    -u  T1-  -^   (l-2n)t^  „/«   N    .        11        n        log(l+Lt)it 

or,  letting  L  =  (x-l)/t  in  the  last  term 

2 

(17)  X  =  1  +  Lt  +  \^^::,^)   H(a^)  +  4(i_^)(x-l)  LI  -  -l^i 

where  H(a,)  is  given  in  equation  (11).  Since  the  logarithmic 
term  in  (l6)  does  not  vary  much,  an  approximate  integration  was 
used  to  obtain  (17). 

Equations  (11),  (15),  (l6),  and  (17)  will  be  used  for 
describing  the  flow  in  region  (J)). 

J>.    Main  Shock  and  Region  (l). 

The  characteristic  equations  (l)  and  (2),  written  with 
pressure  and  velocity  as  dependent  variables,  are 

2 

(18)  p^  +  pau^  +  (u+a)  (p^+pau^)  +  ^^^   =  0 


(19)      Pt  -  pau^  +  (u-a)  (p^-pau^)  +  ^^^   =  0> 


with  p  pressure  and  p  density. 
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These  equations,  when  written  in  this  form,  hold  whether  or 
not  the  flow  is  Isentropic. 

Since  the  positive  characteristics  are  the  ones  coming 
into  the  main  shock,  equation  (l8)  is  used  for  application  of 
Whltham's  rule  (given  in  the  introduction)  for  determination 
of  the  main  shock  motion.   We  write  (l8)  in  differential  fonn 

2 
f^^\  J      J      pa  un  dx  - 

(20)  dp  +  padu  =  -  ii^j^^  —  . 

At  the  shock  the  Hugoniot  relations  are 
2a^  .   _n.      P, 


u  = 


?  (M-M"^);  P  =  -?T  (2tM^-Y+1) 


7+1  ^     ' *   ^        Y+1 


(21) 


_  Po('Y+l)M^   ^   2  _  aQ^(2TM^-Y+l)(M^(T-l)+2) 
(y-1)M^+2   '  (Y+lplI? 


where  M  =  U/a^,  U  is  the  shock  velocity  and  subscript  0 

refers  to  values  ahead  of  the  shock  (region  (O)  Figure  l). 

Substituting  these  relations  into  (20)  we  obtain  a  differential 

relation  between  shock  Mach  nximber  and  its  radius,   x  : 

m 


(22)   -n-^  =  dM|-^^ ^(M^-H) 

L^y^-y-'^        M  v/ (2^^=^-7+1)  (M^(Y-l)+2) 


^m 


+  2M   /(■Y-l)M^+2  ^   M^+1 
M^-1  V  27M^-'Y+1    M(M^-l) 
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This  can  be  integrated  to  give 

(23)   (xj'^(27M^-Y+l)    •  (V(Y-l)(2yM^-Y+l)  +  \/2y(M^(y-1)+2)  j  ^ 


^(y2YM^-Y+l-^('y-l)M^42)^^^p^ 1__  arcsln(^(^'rM^-T+l)-(T-l)(M^(T-l)+2)(] 

N        /^TRT     ]      m^(y+i)^ 


=  constant. 

Using  the  dimensionlesa  coordinates  introduced  in  (3)  we  have  the 
relation  between  t   and  M,   along  the  shock, 

^^r^       dt„  dx^   T  dx^ 

IB 

dx 
Where  — ^  is  defined  in  equation  (22). 

The  motion  of  the  main  shock  as  a  function  of  its  Mach  number 
is  given  by  equations  (23)  and  (24)  in  conjunction  with  em  initial 
condition  specifying  its  position  and  Mach  ntunber  at  a  certain  time. 
Ordinarily,  the  conditions  are  prescribed  at  the  instant  of  deto- 
nation.  If,  for  example,  we  were  considering  the  bursting  of  a 
gas  bubble  at  high  pressure,  the  initial  shock  Mach  number  could 
be  determined  by  use  of  one  dimensional  shock  tube  theory. 
Specifically,  if  it  were  assumed  that  the  gas  inside  the  bubble 
were  steady  with  pressure  and  sound  speed  p^,   and  aj^,   and  that 
the  outside  atmosphere  were  steady  at  p^  and  a^,   the  one 
dimensional  theory  would  give  the  following  relation  which  the 
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initial  shock  Mach  number  must  satisfy. 


Mi 


(25)     ^  ii  .  !^  (M-M-1)  ^    -  (l+^i)M^  -  ^x. 


This  result  Is  derived  in  Appendix  A. 

It  will  be  seen  that  the  only  flow  property  required  to  be 
known  in  region  (1)  is  the  position  of  the  negative  characteristics 
reflected  from  the  main  shock.   These  characteristics  are  defined 
by  the  relation 


dx 


=  u  -  a. 


Since  the  contact  front  follows  very  closely  behind  the  main 
shock,  region  (l)  is  very  small  and  it  suffices  to  approximate 
the  negative  characteristics  there  by  straight  lines.   Effects 
of  the  changing  shock  strength,  due  to  expansion,  are  accounted 
for  by  adjusting  the  slope  of  each  characteristic  to  be  that 
value  of  u  -  a  where  the  characteristic  meets  the  shock.   We 
have,  therefore. 


(26) 


^  =  ^m  +  *'i(M)(t-tJ 


with  w  (m)  =  u  -  a   evaluated  at  the  point  x  ,  t    on  the 
-1-  m   m 

main  shock. 

Equations  (23),  (24)  and  (26)  define  the  main  shock  and 
negative  characteristics  behind  the  shock  as  functions  of  the 
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parameter  M,  which  is  the  shock  Mach  niimber. 

4.  Contact  Front. 

The  contact  front  separates  the  gas  in  region  (2)  from  the 
air  in  region  (1).   The  information  about  the  front  required  for 
the  present  analysis  is  the  location  of  the  front,  and  the  change 
in  slope  of  the  negative  characteristic  upon  crossing  the  front. 
It  will  be  seen  in  the  next  section  that  it  is  necessary  for 
the  determination  of  the  secondary  shock  to  know  the  slope  of 
the  negative  characteristics  in  region  (2).   However,  as  will  be 
shown  below,  in  order  to  determine  this  slope  one  must  know  not 
only  the  characteristic's  slope  in  region  (l)  but  also  the  location 
of  the  point  where  the  characteristic  crosses  the  contact  front. 

In  general,  pressure  and  normal  velocity  are  continuous, 
and  density  and  temperature  may  undergo  Jiimp  discontinuities 
across  contact  fronts.   Since  the  present  analysis  is  concerned 
with  flows  which  are  purely  spherical  or  cylindrical;  i.e.,  with 
one  space  dimension,  contact  fronts  are  always  normal  to  the  flow 
direction,  and  velocity  is  continuous  across  the  front. 

As  the  contact  front  moves  with  the  local  particle  velocity, 
u,   it  will  be  met  by  positive  characteristics,  slope  u  +  a, 
from  region  (2)  and  negative  characteristics,  slope  u  -  a,   from 
region  (1).   In  order  to  determine  the  positive  characteristic 
values  we  shall  assume  equation  (11),  which  gives  — y  a  +  •^  u 
in  region  (5),  to  be  applicable  in  region  (2);  this  due  to  the 
fact  that  the  secondary  shock,  separating  regions  (2)  and  (J>) , 
is  initially  very  weak.   Letting  Q^,  Q^  and  Wg,  w^^  denote 
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_1^  a+iu  and  u-a  In  regions  (2)  and  (l)  respectively, 
7-1      2 

We  have,  due  to  continuous  velocity  at  the  contact  front 


Q2  -^  Y^  *'2   =  ^1  -^  T^  ^1 
(27)         <(    or 

W2  -  w^  =  -  (t-i)(V'^i^  • 

For  Isentroplc  flows  we  have  the  relation  between  pressure  and 
sound  speed  [Ref o5jP«95] * 

Pi   ^a^^ 

where  subscript  1  stands  for  some  initial  state.   Consequently, 
due  to  continuous  pressure  at  the  contact  front 


1       1 


(28)         y  or 


^12 


2^2  -   «2  =  i:^(2Q^.w^) 


(It  should  be  noted  that  in  the  derivation  of  equations  (27)  and 
(28)  the  assumption  has  been  made  that  7,   the  ratio  of  specific 
heats,  has  the  same  magnitude  in  both  regions  (l)  and  (2).   The 
extension  to  different  7's  is  direct  but  not  too  simple,)   We 
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obtain  the  slope  of  the  negative  characteristics  in  region  (2) 
in  terms  of  kn^ 
(27)  and  (28). 


in  terms  of  known  quantities  by  eliminating  Q^   from  equations 


(29)      w, 


Hiiii    +  2(1  -  ^)Q. 

^  ,  .   2  ^ 

r-1  a,^ 


Evaluation  of  (29)  is  not  direct,  since  w,   is  knovm  only 
as  a  function  of  main  shock  Mach  nvunber  (eqn.(26))  and  Q^  is 
knovm  only  as  a  point  function  of  x  and  t   (eqn.(ll).   We 
must  determine,  therefore,  the   (x,t)   coordinates  of  the  contact 
front  in  terms  of  the  parameter  M. 

First  we  express  a,   in  terms  of  w,   and  Qp  : 
velocity  is  continuous  across  the  contact  front 


since 


•^  -  '^i  =  rT  ^2 


+  a. 


combining  this  with  the  first  of  equations  (28) 

!Q2  -  «! 


(50)  a^  = 


1  +_2_!21 
7-1  a,. 


Next,  let  X  =  c(t  )  be  the  equation  of  the  contact  front. 
The  negative  characteristic  given  in  equation  (26)  meets  the 
contact  front  at  the  point   (x  , t  ).   Hence,  at  this  point 
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c(tc)  =^m  +  ^(''Jt^c  -  ^m)- 


Differentiating  with  respect  to  t   we  have 


dc 
^0 


'"i^H^^f'i^-itV*™.)  -  Vij- 


Primes  indicate  differentiation  with  respect  to  M.   As  the 
contact  front  moves  with  local  particle  velocity. 


dc    ^^c 
c     c 


Combining  the  two  above  equations,  using  w  =  u  -  a 

(51)      a,  = -^  |x'  +  w'(t  =t„)  -  w^t'  f   . 
^^   '  1   dt   ?  m    1^  c  m'    1  m  j 


Eliminating  a^  from  (50)  and  (51)  we  obtain  a  differential 
equation  relating  t   and  M 

dt    f^-^^]   K  -^  n^S-^m)  -  ^^; 
(32)        ^   ^        ^^' 


^T»r 


20^  -  w^ 


Equation  (32)  together  with 

(55)  ^c  =^m  +  ^(M)(t^-tJ 
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define  the  coordinates  of  the  contact  front   (x  ,t  )   as 

*  c'  c' 

functions  of  the  parameter  M.   This  result  enables  us  to 
evaluate  the  change  In  slope  of  the  negative  characteristic 
upon  crossing  the  contact  front  as  given  In  (29).   In  fact, 
we  have  directly  the  negative  characteristics  In  region  (2) 

(54)  ^  =  x^  +  (t-t^)w2  , 

where  t   and  x   are  given  by  (32)  and  (35)*  and  w„  Is 
given  by  (29). 

5.  Secondary  Shock. 

We  now  have  enough  Information  about  the  flow  to  fit  the 
secondary  shock  between  regions  (3)  and  (2).   It  Is  seen  that 
those  negative  characteristics  arising  In  the  expansion  region 
(3)  tend  to  fan  Into  the  Increasing  x  direction.   Those  which 
arise  at  the  main  shock  will  tend  to  point  more  and  more  toward 
the  decreasing  x  direction.   This  Is  due  to  the  fact  that  the 
shock  weakens  as  It  expands,  as  can  be  Inferred  from  equations 
(22)  or  (23).   Thus,  as  the  shock  weakens  the  quantity  u^^  -  a^, 
which  is  the  characteristic  slope,  decreases. 

The  secondary  shock  must  be  Inserted  so  that  the  shock 
relations  are  satisfied.   For  a  weak  shock  the  shock  conditions 
show  that  Its  slope  at  each  point  In  the   (x,t)  plane  Is  the 
average  of  the  slopes  of  the  characteristics  meeting  It.   Since 
we  have  the  equations,  (17)  and  (34),  describing  the  characteristics 
which  meet  the  shock,  we  can,  using  a  technique  developed  by 
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Whltham  [Ref.9],  obtain  the  shock  motion. 

The  characteristics  which  arise  In  the  expansion  fan 
and  at  the  main  shock  are  respectively 

X  =  1  +  Lt  +  f (x,t) 
X  =  x^  +  WgCt-t^) 

where  f(x,t)   represents  terms  appearing  In  equation  (17)» 

and  X  ,  t  ,   and  w^  are  fvinctlons  of  main  shock  Mach  number, 
c    c         d 

M,   as  given  In  equations  (33) »  (32)  and  (29).   Assxomlng  the 
secondary  shock  path  to  be  given  by 


Xg  =  1  +  S(t3) 


we  have  along  the  shock  path 

(56)  S(tg)  =  W2(tg-t^)  +  x^  -  1 

=  Ltg  -f  f  (x^+(t3-t^)w2,tg)   . 

Expressing  the  shock  slope  as  the  average  of  the  characteristic 
slopes  as  obtained  from  (35), 


dS    1  ,     ^  "^  ^t, 


^^  =  ^  (w^  +  j-T-f-)  • 


Letting  the  parameters  M  and  L  and  the  variable  x   all  be 

s 
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dependent  on  t_,  we  can  obtain  another  representation  for 

^|-  ,  using  (35)  and  (36) j  letting  primes  represent  differentiation 
with  respect  to  M: 


0  o 


This  equation^,  when  combined  with  the  previous  expression,  yields 


(37)      ^m--    (1  "  ^x)(^i  +  '^^(^s-tc)  -  ^a^c) 


Eliminating  L  from  (37)  by  use  of  (36)   we  obtain,  after  some 
algebra,  a  relation  between  t„   and  M: 


(38)     ^   (t^  -  1  °  ^c-'a'  ).[,,-!-  Vgl-(-  2W.  .  1^  (f) 

s  s 

+  tf  K  +*'2(V^c)  -  Vc)^ 


Here,  f  and  f   are  evaluated  at  x  =  x  +  (tg-t  )w2,  t  =  t^. 
Integrating  both  sides  of  (38)  with  respect  to  M  we  have  finally 

(39)   [x^  -  1  -  t^Wg]^  =  tg  ^"[x^  -  1  -  t^W2].[-  2w.  +  y^) 

f. 


+  t;  (^c  +  ^2(V^)  -  ^2^c)3^« 
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where  M   Is  the  initial  Mach  nvunber  of  the  main  shock  wave. 
Equation  (59)  is  rather  messy  and  must  be  solved  by  an 
iterative  numerical  technique.   Since  for  a  fixed  M  the 
right  hand  side  of  (59)  increases  raonotonically  with  t^  and 
the  left  hand  side  is  a  function  of   M  alone,  it  is  possible, 
for  each  increment  in  M,   to  solve  for  the  corresponding  t^ 
by  a  simple  "trial  and  error"  scheme.   Once  t^  is  determined, 
x„  is  given  by 

(40)  ^s  =  ^c  +  *2(V^c)' 

The  technique  described  above,  while  useful  to  describe 
the  initial  motion  of  the  secondary  shock,  is  valid  only  when 
this  shock  is  still  weak.   In  the  present  situation,  however,  this 
shock  strengthens  as  it  is  carried  outward  by  the  expanding  gases 
and  another  method  must  be  used  to  describe  its  subsequent  motion. 
(The  change-over  point  for  a  specific  problem  was  arbitrarily 
chosen  to  be  either  when  the  shock  begins  to  turn  back  toward 
the  origin  or  where  the  shock  strength  becomes  1  in  magnitude, 
whichever  comes  first.   Shock  strength  is  defined  in  this  case 
as  Mg  -  1,   where  M^  is  the  shock  Jfech  number.) 

For  moderate  to  strong  shocks,  the  technique  derived  in 
Section  3  for  describing  the  main  blast  wave  is  applicable. 
It  must  be  altered  in  the  present  case,  since  the  flow  coming 
into  the  secondary  shock  is  non-uniform  and  terms  describing  this 
incoming  flow  variation  along  the  shock  must  be  Included.   To  do 
this  we  start  with  the  characteristic  differential  relation 
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corresponding  to  (20),  but  this  time,  since  the  negative 
characteristics  are  the  ones  coming  into  the  shock,  we  use 


(^1)  dp  -  padu  =-£^^  . 


As  the  expansion  region  (3)  is  one  of  isentropic  flow. 

The  Hugoniot  relations  given  in  (21)  are  used  with  the  variables 
a^,  Pq,  Pq  replaced  by  a,,  p,,  p,.   Also  the  velocity  relation 
is  altered  to  account  for  the  fact  that  the  shock  is  backward 
facing  and  the  flow  in  region  (3)  is  non-stationary,  giving 

^  =  ^5  -  ^  (M-M  ^) 

where  M  «  (u-,-U  )/a^,   and  U   is  the  secondary  shock  velocity. 
Substituting  these  relations  into  (4l)  we  obtain  equations  which 
hold  along  the  secondary  shock 
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{fi(^^24p)-(4i^^)< 


,m2-1  .   1  -^i^W^^^s 

a  I 


+  ( 


2a  (M^-l)-('Y+l)Mu 


(Y+l)Ma,-2a^(M^«l)-R 


1^3-  ^  - 


^) 


?5^ 


2M      M^+l 
MR 


2'yM  -7+1 


ax    s 


U  =  u,  -  Ma, 
s    3     3 


with  R  =    \/(2'YM^-'y+l)(l^(7-l)+2) 


The  quantities  u^,  C,,   and  their  derivatives  are  obtained  from 
equations  (11)  and  (l6) .   Equations  (42)  are  solved  stepwise  from 
a  given  point  where  x,  t,   and  U„  or  M  are  prescribed. 


6.  Application  to  a  Specific  Problem,, 

We  shall  describe  here  the  application  of  the  present  theory 
to  a  specific  problem  which  has  been  solved  theoretically  by  Erode 
[Ref.2],  and  carried  out  experimentally  by  Boyer  [Ref.l].   A 
sphere  of  radius  1,  containing  a  gas  at  22  atm.,  and  at  a  temperature 
of  299  K,  is  assumed  to  be  surrounded  by  air  at  1  atm.  and  the  same 
temperature.   The  compressed  gas  is  also  air,  and  the  ratio  of 
specific  heats,   y,   is  assumed  to  equal  1.4  everywhere.   The 
initial  Mach  number  of  the  main  shock,  as  determined  by  one 
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dimensional  shock  tube  theory  (Appendix  A),  Is  1.846.   With  this, 
the  constant  on  the  right  hand  side  of  the  main  shock  equation 
(25)  Is  26.1.   For  this  case  the  sound  speed  and  particle  velocity, 
made  dlmenslonless  with  respect  to  the  flow  In  region  (O),  are  as 
follows,  at  the  Initial  point  x  =  1,  t  =  0: 


Uq  =  0 

^0  =  1 

Uj^   =  1.087                 ^ 

a^^   =  1.252 

Ug   =1.087 

ag   =      .729 

"5 '  id  *  ^) 

a,  =  ^(5  -  ^) 

U4    =   0 

a^  =  1. 

The  leading  characteristic  of  the  expansion  fan  In  region  (5)  Is 
(x-l)/t  =  -1. 

For  region  (3)  we  have,  corresponding  to  equations  (11), 
(16)  and  (17) 

(6.11)    2.5aj  +  .5u^  =  2.5  -  ^^(5  -  ^)  (1  +  ^f 

(6.16)    u^  -  a^  =  -^  . [(1  +  -^)  -  3^(1  -  -^^)] 


(6.17)    X  .  1  +  Lt  -  ""''^'"'^'"^'"'[(1  .  i^)  -  i§r(l  -  ^^)] 


The  contact  front  Is 
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5a^  +  u^  -  w^ 


(6.35)         ^c  =^m  +  *'i(W)(VtJ, 

w,  (M)   is  the  negative  characteristic  slope,   u  -  a,   in  region 
(1),  It  Is  evaluated,  as  a  function  of  the  main  shock  Mach 
number  M,   by  using  the  Hugonlot  relations  In  equation  (21). 
In  region  (2)  the  negative  characteristic  slope,   Wg,   is 

(6.29)         Wg  =  .909w^  +  .182(2. 5a^  +  .5u^) 

In  Figure  2  the  present  theory  Is  compared  with  that  of 
Boyer  and  Erode.   The  results  of  several  experimental  explosions 
carried  out  by  Boyer  are  indicated  by  dashed  lines.   Erode 's 
work,  numerically  integrating  the  Lagranglan  equations  of  motion, 
is  indicated  by  solid  lines.   The  present  theory  is  given  by  the 
curves  composed  of  heavy  dots.   A  simplified  determination  of  the 
secondary  shock  is  Indicated  by  x's.   For  this  the  one-dimensional 
solution  was  used  in  region  (3)-   This  is  given  by  omitting  the 
last  term  in  equations  (6.11),  (6.l6)  and  (6.17),  and  also  setting 
f  and  f   equal  to  zero  in  the  secondary  shock  equation  (39). 
The  time  to  compute  the  two  shocks  and  the  contact  front  on  an 
IBM  704  was  about  3  minutes. 

Since  Boyer 's  experimental  work  involved  the  bursting  of 
a  glass  sphere  the  differences  between  his  result  and  Erode 's 
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may  be  due  to  the  presence  of  glass  particles  in  the  flow  field. 
As  the  present  model,  a  hypothetical  sudden  expansion  of  a  gas 
sphere,  is  the  same  as  Erode' s  it  is  expected  that  the  results 
be  comparable  with  his.   The  agreement  as  to  the  initial  motion 
of  the  main  shock  and  contact  front  is  quite  good.   Brode  shows 
the  secondary  shock  expanding  slightly  more  and  imploding  sooner 
than  both  the  present  theory  and  Boyer's  experiments.   The 
differences  in  the  initial  motion  could  be  due  to  the  fact  that 
the  shock  is  initially  weak.   Brode,  using  the  von  Neumann- Rlchtmyer 
pseudo  viscosity  technique,  which  spreads  out  the  shock,  might  have 
difficulty  in  determining  the  true  location  of  a  weak  shock.   Boyer's 
experimental  technique  which  used  the  schlieren  shadowgraph  to 
determine  the  secondary  shock,  has  difficulty  in  locating  weak 
shocks  which  give  correspondingly  weak  shadows.   As  to  the  determi- 
nation of  the  time  of  implosion  Erode 's  use  of  Lagranglan  mass 
coordinates  may  lead  to  inaccuracies  near  the  flow  center  as  the 
density  there  is  quite  low.   Most  of  the  mass  is  located  near 
the  contact  front  and  main  shock.   However,  since  the  present 
theory  becomes  less  valid  as  the  distance  to  the  initial  point, 
X  =  1,  t  =0,   increases}  it  is  rather  difficult  to  precisely 
determine  the  causes  of  the  small  discrepancies  in  the  different 
results. 

The  approximation  using  the  one-dimensional  flow  solution 
to  determine  the  secondary  shock  agrees  remarkably  well  with  Boyer's 
experimental  work.   This  approximation  may  be  useful  as  it  gives  a 
qualitative  description  of  the  flow  and  affords  a  marked  simplification 
of  some  of  the  equations.   Overall,  however,  the  agreement  between  the 

different  theories  is  quite  good. 
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Appendix  A.  One  Dimensional  Shock  Tube  Theory. 

Consider  a  tube  separated  Into  two  regions  (Figure  5) >  each 
region  containing  a  uniform  stationary  gas  but  at  different 
pressures.   The  separation  point  Is  assiimed  to  be  at  x  =»  0,  of 
a  coordinate  along  the  tube,  and  the  high  pressure  region  Is  In 
X  <  Oo   At  time  t  =»  0  the  separation  Is  removed  and  a  shock 
moves  Into  the  uniform  region  (0).   Regions  (1)  and  (2),  separated 
by  a  contact  front,  are  both  iinlform  and  have  the  same  pressure 
and  particle  velocity.   Region  (5)  Is  a  centered  expansion 
separating  the  uniform  regions  (4)  and  (2).   Knowing  the  flow 
properties  In  regions  (0)  and  (4)  It  Is  desired  to  determine 
the  shock  strength  and  the  slope  of  the  last  characteristic  of 
region  (5) J  I.e.,  the  characteristic  separating  regions  (2)  and  (3). 

In  the  expansion  region  the  solution  Is  [Ref .5*p.l04] 


(A.l) 


,  =  .(i^  a,  »  f) 

^  -  (l-n)(ai^  +  |) 

i-ji 

where  p.  =  ^y  .   The  Hugonlot  shock  relations  giving  u,   and 


p,   are 


(A.2)     u^  =  (l-n)aQ(M-M''^),  p^  =  Pq((1+h)m2-h) 


Since  the  velocity  and  pressure  are  continuous  across  the  contact 
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frontj,  Up  and  p^  are  also  given  by  (A. 2).  Letting  S 

denote  the  slope  of  the  last  characteristic  of  the  expansion 

region  (5),  we  have  at  this  characteristic  due  to  continuous 
velocity  and  pressure 


-1, 


a^^  +  S  =  aQ(M-M~  ) 


(A. 3)  J 


V 


P4 


l+M. 
S))  ^ 


P^((1+H)*P-M.) 


Eliminating  S  in  (A. 3)  an  equation  for  M  is  obtained 
(A.4)     ^(1  -  ^(M-M-^))  ^     =  {l+n)M2-^.. 


This  gives  the  required  shock  Mach  number  in  terms  of  the  initial 
conditions. 
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Figure  2.   Experimental  and  Theoretical 
Spherical  Blast  Results 
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